Bat dynamics modelling as a tool for conservation management in subterranean environments

Bat species inhabit subterranean environments (e.g., caves and mines) in small areas with specific microclimatic conditions, during various periods of their life cycle. Bats can be negatively influenced by microclimatic changes within their roosts if optimal habitat patches become unavailable. Therefore, proper management solutions must be applied for the conservation of vulnerable bat populations, especially in show caves. We have pursued an ensemble species distribution modelling approach in subterranean environments to identify sensible patches for bats. Using multi-annual temperature monitoring and bat distribution surveys performed within ten caves and mines, including show caves, we modelled relevant habitat patches for five bat species. The temperature-based variables generated from this approach proved to be effective when processed via species distribution models, which generated optimal validation results, even for bats that were heavily clustered in colonies. Management measures are proposed for each show cave to help in long-time conservation of hibernation and maternity colonies. These measures include creating suitable microclimatic patches within the caves by ecological reconstruction measures, tourist management practices in relation to bats, and show cave fitting recommendations. This approach has never been performed at this scale due to the complex geostatistical challenges involving subterranean environment mapping and can be further used as best practice guidelines for future conservation projects.


Introduction
Subterranean environments (SE: caves, mines) are one of the most widespread habitats worldwide and harbour diverse ecosystems with a high degree of species endemism [1]. Although a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 During a hibernation season, they can exit torpor and move several times within the roost or between a cluster of roosts in search of favourable climatic conditions [33]. The animals can increase their arousal episodes and even forage near the roost if the air temperature reaches a certain threshold (e.g., approximately 10˚C for R. ferrumequinum [27]).
Anthropogenic factors can cause disturbances for hibernating or nursing bats, often expressed as mass tourism. These interventions can lead to strong disturbances in the energy equilibrium of torpid bats, such as weight loss [34], which can lower the survival rate of the exposed population [35] or the abandonment of optimal nursing sites [36]. Open space cavedwelling bats are especially vulnerable to external factors which can cause arousals, such as artificial lighting, changes in air circulation, and slight increases in air temperature, which mainly occur during visiting hours [37,38]. On a global scale, climatic changes may also influence air circulation in caves, leading to increased air temperatures and loss of suitable habitats for a broad group of cave-dwelling species [39]. The recent climatic evolution correlated with anthropogenic factors amplifies the need for extensive spatial biodiversity and microclimatic monitoring within SEs, as a framework for wildlife and habitat conservation [40].
Desipite the vast distribution of SEs worldwide and their large species diversity [1], spatial distribution studies regarding subteranean habitats remain scarce. Species distribution models (SDM) can help predict potentially suitable habitats within a site and can be used as conservation and management tools [41], but SDMs have rarely been applied to SEs due to a series of challenges in building environmental variables [42]. The existing studies presented limited results for the spatial microclimatic preferences of animals in SEs, often for a single cave or for a short time frame [43,44]. Some approaches focused on collecting occurrences from within the sites and projecting them on external environmental variables [41,45,46], allowing for an increased number of analysed sites at the expense of detailed habitat mapping for each particular SE. Cave-dwelling bats use specific areas within a SE for shelter or reproduction [47] and are conditioned by a limited set of environmental factors, mostly revolving around stable microclimates.
Although bats are one of the most commonly studied taxons in SEs, especially within the temperate regions [48], conservation efforts seem to have mixed results, pointing out that more systematic site specific studies are required, or that more complex disitrbution recognition patterns of their roosting preferences need to be undertaken to implement successful conservation measures.
Given these limitations, our goal was to map the distribution of several cave-dwelling bat species during multiple key periods of their life cycle and to understand how the animals use the sites to propose optimal management solutions, especially in show caves. All bat species within our studied range are strictly protected and most of the selected roosts are included within the EUROBATS important underground sites list. The species are included in the European Natura 2000 Network, and therefore require specific conservation measures to ensure population connectivity and habitat integrity, including the creation of new protected natural areas [49].
We used species distribution models (SDM), specifically ensemble SDMs (ESDM) to reduce biases of individual modelling techniques by averaging predictions across a multitude of models [50]. It is generally accepted that this approach provides a better predictive performance compared to a single model running algorithm [51][52][53]. We created a series of temperaturebased variables extracted from a systematic air temperature monitoring effort. We further compared the results to the observed anthropogenic impact for each site to identify relevant conservation and management measures and help reduce the impact on bat populations. We hypothesized that bats' spatial and temporal distribution within SEs can be accurately predicted via SDMs, aiding future management efforts. This novel approach was challenged by a series of limiting factors associated with creating subterranean environmental spatial datasets and reliable occurrence sampling strategies. Results can be used as a conservation tool for bat colonies that inhabit caves and mines in most bioregion, especially show caves, aiding managers in their conservation efforts.

Methods
This study was performed in compliance with the recommendations described within the Bat Surveys Good Practice Guidelines of the Bat Conservation Trust. The field protocol was authorised by the Romanian Academy-Natural Monuments Commission (#3660/22.11.2012). The data collection efforts were designed to reduce bat arrousals in all the monitoring periods, with a reduced number of field surveyors per site and a short observation time frame near the animals.

Data collection
The study was performed in Romania, a country with diverse karst landscapes and large cavedwelling bat populations [54]. We chose ten SEs (Fig 1, Table 1) located in the steppe (n = 3), continental (n = 1), and alpine bioregions (n = 6). We classified them into three categories, based on their origin and current use: wild caves, show caves, and mines.
The number of bats in colonies was determined using flash photography to minimize temperature-based arousals [56]. Observations were limited to the extent of the show cave paths, where most bats roosted (SEs 4, 5, 6, and 10).
Data collection was split into three periods: hibernation start (SO: September-October), mid to end hibernation (NM: November-March), maternity start (AM: April-May). Each period was sampled three times per year (beginning, mid, end), and each SE was observed for two consecutive years, between 2010 and 2016 ( Table 1). The 2014-2016 period did not cover the complete maternity interval to minimize human impact on bat populations. The surveys from this period were useful for model development, as bats already form the colonies at the start of the maternity period.
The spatial extent of the studied SEs, which represented the modelling boundary, was modified from published maps (Fig 2), but we also performed new surveys where the data was outdated or inexistent (SE 1-3). We used a laser meter (Leica DistoX) fitted with an inclinometer and PDA-Paperless Cave Survey software [57,58] for all cave surveys. The SEs which had published maps were also surveyed, but only using a central line which helped us georeference the existing maps. The survey data was processed using Compass Cave Survey and Mapping Software. We merged the final maps into a single layout, but each SE had its own spatial reference. This was used only as a visual aid for results presentation, while the models were performed at their original scale of all SEs.
The focal bat species were Rhinolophus ferrumequinum, R, hipposideros, R. euryale, Myotis myotis / M. blythii, and Miniopterus schreibersii. It was impossible to distinguish between colonies of M. myotis and M. blythii without causing arousal, and therefore they were treated as a single group. Less frequent species, such as R. mehelyi, Nyctalus noctula, and Pipistrellus pipsitrellus (Table 2), were also recorded.
Mapping bat occurrences and human disturbances within the SEs were performed using the same survey method. We created a central survey line, discretely marking stations for future reference. We then used these stations to determine the location of the occurrences. We recorded air temperatures close to the animals each time we collected occurrences. Supplementary data include the height at which bats were located relative to the floor and the distance from the nearest entrance (spot variables). The distance was extracted via Network Analyst extension (ESRI). Human activity within show caves (SEs 4, 5, 6, 10) and a closed-circuit cave (SE 9) was measured with the help of people counters using a method described in [64,65]. We used the number of passes recorded by people counters (tourist passes) as a proxy for the number of tourists who visited the caves. This proxy does not reflect the number of tourists because people counters log each movement regardless of direction. Non-touristic SEs were not monitored using this method, but we attributed 100 human passes per year for scalability. The resulting spatial datasets were processed within the ArcGIS (ESRI) environment.

Air temperature monitoring and environmental datasets
We measured hourly air temperature values using two methods: continuous monitoring via data loggers and spot temperature measurements with instant probes, as described in [65]. The time frame of the temperature monitoring study overlapped the species observation intervals. The climatic monitoring points were positioned in the field using the above survey method. The continuous temperature measurements were performed using Gemini Tinytag Plus 2 loggers (±0.01˚C accuracy, 0.02˚C resolution) and Hobo Pendant Temperature loggers (±0.53˚C accuracy, 0.14˚C resolution). Before the release of these data loggers, we used iButton Hygrocron from Maxim Integrated (±0.5˚C accuracy, 0.0625˚C resolution with post-processing corrections) for half of the sites (SEs 1, 2, 3, 7, 8). The loggers were positioned 2 m above the cave floor, in the central area of the passage section or chamber. Spot air temperature measurements were collected via a Vaisala HMP70 temperature probe (±0.2˚C accuracy, 0.01˚C resolution), both close to isolated bats or colonies and at specific points where site morphology could influence air temperature values. Spot measurements were later used to extract statistical data regarding bat preferences for an independent validation dataset in the temperature analysis. External air temperature recordings were collected near the SEs entrances, using meteorological stations (SE 5,6,9) or data loggers.
Cave air temperature surfaces were interpolated (0.5 m resolution) using the data logger information (S1 Table in S1 File) via the natural neighbour method [66]. Cave walls were used as boundaries. We used a non-parametric interpolation, opposed to the geostatistical approaches previously used, such as kriging [43], because the geometry of most SEs does not allow for a gridded sampling strategy. Daily interpolations were created, and then averaged to capture temperature variations for each relevant bat activity period. The observations were performed in the middle of the galleries; however, to interpolate a valid surface, each point was copied two times, placing the additional points perpendicular to the cave walls outside their boundaries. Given the fact that Meziad show cave (SE 10) stretches on two main levels (LVL1 and LVL2), which slightly overlap but have different microclimatic regimes, we analysed them as two distinct SEs. The interpolated results were merged into one layer, prioritizing the lower level because it held most of the collected occurrences.
Two validation methods (dependent and independent) were used to check the errors of the interpolations. This was achieved using cross-validation with the interpolation dataset, created   [67]. We included data on human-induced temperature changes (resulting from tourist presence and light systems) collected from show caves used in our study [65].

Statistical analysis
To assess whether there are statistically significant differences between the total number of bats from all the SE over the three analysed periods, the number of individuals from each species during the activity periods, and the total number of bats per each SE over the three activity periods, we performed a series of Friedman's tests. The effect size of the test was measured with Kendall's W coefficient. Pairwise comparisons using paired Wilcoxon signed-rank tests were used to identify the activity periods with significant differences. We, furthermore, calculated the Spearman rank correlation coefficient to estimate the measure of association between the cave length and size of bat populations and also between the cave length and species diversity during the analysed activity periods.
We used Generalised Linear Models (GLM) with stepwise backward elimination to evaluate the differences between five bat species (Rhinolophus euryale, R. hipposideros, R. ferrumequinum, Miniopterus schreibersii, and Myotis myotis / blythii) in relation to torpor height relative to the SE floor, the distance from the entrance, and measured air temperature for all the monitoring periods. We performed a GLM with a binomial distribution using predictor variables for each species. We tested the variable multicollinearity using the Variance Inflation Factors (VIF) function in R, eliminating values greater than 10 [68]. The best model selection was achieved using the lowest Akaike's Information Criterion (AIC) value. To identify the relation between bat abundance per species and periods, compared to the number of people who visited the show caves, and for comparisons of cave metrics and bat populations, we performed Spearman's correlation tests. The analyses were performed in R statistical software version 3.6.3 with the "stats" package [69].

Species distribution modelling
We classified occurrences per species into three periods, as mentioned above, and excluded species with less than 25 occurrences [70]. Spatial autocorrelation revealed heavily clustered datasets (Moran's I index = 0.72, z = 926.14, P � 0.001) because bats form clustered colonies to conserve energy [24]. To account for the autocorrelation errors caused by the species' clustering during torpor or maternity, we have chosen to model the distribution of each species in all SEs for each activity period with no spatial thinning, increasing model variability. Some bat species were not found in all of the analysed SEs (Table 2). Our models were calibrated for each species and selected period within the geometry and environmental datasets of all SEs.
We used ESDM (Ensemble Species Distribution Models) from the SSDM (Stacked Species Distribution Models) R package [71], generating an ensemble model from multiple algorithms: classification tree analysis, multivariate adaptive regression splines, generalized linear models, artificial neural networks, general additive models, support vector machines, and Maxent. The ensemble method generates robust outputs using a large set of algorithms, improving the accuracy of the predictions [51]. These were stacked together into one AUC weighted projection. Models were calibrated with a random subset extracted from the occurrence data (75%) and evaluated with the remaining 25% test datasets [72,73]. The variable contribution was calculated using the Pearson correlation coefficient, SSDM R package. The performance of the models was evaluated using the Receiver Operating Characteristic (ROC) Area Under the Curve (AUC) [74] and Cohen's K test [75]. Both range from 0 to 1, and the maximum value indicates the best fit of the model. In addition, we used a ten-percentile omission error, obtained with the method described by Gherghel and Martin [76], because it is an independent metric extracted from the validation dataset and shows how well the model discriminates between presence and absence predictions. We summed the final predictions using ArcGIS (ESRI) cell statistics into separate datasets per species and period and further summed these into a cumulative dataset for each SE with no prioritization, considering that all bat species are equally sensitive to human disturbance.

Data collection
SEs varied from large wild caves with complex climatic systems (SEs 7, 10; Table 1) to small horizontal artificial sites (SEs 2, 3) with strong seasonal climatic shifts.
Bat distribution compared to the SE temperature range per activity period was plotted in the supplementary materials (S1-S3 Figs in S1 File). Most bats were found in more stable temperature ranges during hibernation and in more variable temperature ranges during the maternity period.
The five focal bat species were recorded in most of the analysed SEs, while the less frequent bat species were observed only in a few sites (Table 2).

Air temperature monitoring and environmental datasets
Most of the SEs presented low and tolerable errors for bat activity (MAE < 0.4, RMSE < 0.4), except for Gura Dobrogei and Valea Cetăţii caves (SEs 1 and 4) which showed higher RMSE values (0.6 and 0.9, respectively). These were recorded in subterranean sectors which were located closer or between the entrances, but also where longitudinal profiles show vertical changes. Slight differences in accuracy and resolution of the equipment also contributed to some errors. The values were recorded in patches with low bat abundances and were insignificant compared to the dynamic variations recorded at the entrance of the SEs or the ecological requirements of bats, therefore, we considered the errors negligible. The independent validation also showed minor differences between the measured and the predicted values, with a maximum MAE of 0.4 and RMSE of 0.9 for Gura Dobrogei cave (SE 1), validating the interpolations.
Subterranean environments have a vertical temperature gradient, with hotter air concentrating closer to the ceiling. Due to field collection restrictions, this was not covered by the data logger monitoring method, especially in high chambers. Nevertheless, some spot measurements were collected at much greater heights, matching the height of the bats. Differences are shown in the independent test datasets (S2 Table in S1 File), but they do not exceed the tolerable errors relevant for torpid or nursing bats (< 0.9˚C).
Some temperature disturbances were recorded within the show caves, with an increase by 0.4˚C near the lights in Polovragi cave (SE6). These appeared during the visiting period and were proportional to the number of people recorded on the people counters. They reverted to the initial stable temperature (IST) in approximately 40 hours. Muierilor cave (SE 5) recorded higher temperature spikes in the touristic season (2˚C, over 250000 passes), which slowly reverted to IST in a month. Increased values of 0.5 to 1˚C were still recorded in medium-sized halls during the winter, where smaller groups of tourists were stationed. During winter, the number of tourists dropped; therefore, temperature values reverted faster to the IST (one day).
Valea Cetăţii and Meziad caves (SEs 4 and 10) showed minimal temperature variations associated with human disturbances.

Statistical analysis based on spot measurements and observations
Friedman's test showed that the number of bats was significantly different during the three analyzed periods (P < 0.05), with a moderate effect size (W = 0.44). The results of the Wilcoxon signed-rank test showed significant pairwise differences between NM-SO (P < 0.05) and NM-AM (P < 0.05). While testing for differences between the number of individuals from each species and activity periods, the results of Friedman's test showed significance for Myotis myotis/Myotis blythii (P = 0.0368), Rhinolophus ferrumequinum (P = 0.0003), and R. hipposideros (P = 0.0129). The following SEs registered significant differences among the number of bats over the three activity periods: Muierilor cave (SE 5: P = 0.0498, W = 0.75), Polovragi cave (SE 6: P = 0.0224, W = 0.95), and Canara Tunnel (SE 2: P = 0.0067, W = 1). At the beginning of the hibernation period (SO), the cave length was strongly correlated with the size of the bat populations (Spearman's rho = 0.9192, P < 0.005), but this correlation did not follow during the NM and AM periods. Species diversity also correlated positively with cave length during SO (Spearman's rho = 0.84966, P < 0.005), but not during NM and AM.
Results of the GLM showed that the species were significantly influenced by height, distance from the cave entrance, and measured temperature (spot variables-S3 Table in S1 File) during at least one of the analysed activity periods (Table 3).
At the beginning of hibernation (SO), R. hipposideros, R. ferrumequinum, M. schreibersii, and M. myotis / blythii showed a preference for the higher temperatures, closer distances from the entrance (except for M. myotis / blythii), and higher hibernation heights (only R. euryale and M. schreibersii). In contrast, during the mid-end part of the hibernation (NM), the species preferred lower heights, greater distances from the entrances, and lower temperatures. During maternity (AM), all species except R. ferrumequinum searched for higher temperatures, while R. ferrumequinum and M. myotis / blythii also showed a preference for greater heights relative to the SE floor (Table 3).
Correlation between bat abundance per species and observation periods compared to the number of tourists' passes in show caves showed insignificant values, except for R. hipposideros during the start of the hibernation (Table 4).
Thermophilic species such as R. mehelyi and R. euryale searched for warmer areas during hibernation compared to crevice-dwelling species, such as Nyctalus sp. or Pipistrellus sp. Most Rhinolophus spp. choose roosts with lower temperature variations and preferably one entrance. Rhinolophus spp. hibernated for more extended periods than crevice-dwelling species, which preferred cracks in large chambers (in SE 10), with higher temperature amplitude. Distance from the entrances increased when exterior temperature decreased, except for M. myotis / blythii, which got closer to the entrance in those conditions. Most species decreased their roosting height towards the end of the hibernation in search of lower temperatures to conserve energy. In contrast, greater heights were linked to warmer areas optimal for nursing during maternity for all the studied species.

Species distribution models
Results showed potential suitable habitats for most of the focal species in most SEs, but not all species have yet to colonize those areas. Ensemble models performed well, with a minimum of 0.75 and an average of 0.89 training AUC. Cohen's Kappa also showed optimal averaged results, reaching 0.81. Suitable habitats per species were accurately predicted, with an averaged 10% omission rate of 8.3%. The selected model variables changed for each activity period according to each species' environmental requirements. During the start of hibernation, minimum temperatures were more relevant; during the mid-end hibernation, the maximum temperatures and the temperature ranges were significant, while during the maternity period, the minimum and maximum temperatures were mainly considered ( Table 5).
The sum of the binary models per period (Figs 3-5) showed where an aggregation of species was more likely to appear. During SO, the models showed a higher abundance of optimal habitats in areas that mostly overlapped the static temperature sectors of the SEs. Colder SEs presented higher SDM values, such as SE 4, SE 7, and SE 8. These favourable sectors became more restrictive during NM and concentrated at lower stable temperatures than SO, but with some degree of temperature variability. The colder SEs lost some degree of favourability compared to the more complex roosts, such as SE 5, SE 9, and SE 10, where the animals could find optimal climates for torpor. The AM period showed substantial bat abundance and distribution changes, with essential sectors closer to the entrance. Meziad cave (SE 10) showed a climatic inversion, where hot air reached and accumulated in the upper level, at the far end of the show cave sector. Nevertheless, the SDMs accurately predicted these important roosting areas for bats (M. myotis). Cold caves such as Valea Cetăţii or Stogu (SE 4, SE7, respectively) had little to no predicted SDMs during AM.

Discussions
This study emphasized the importance of using ensemble species distribution models on subterranean environments to predict suitable habitats for bat species during their lifecycle. For this novel method in terms of scale and functionality, we applied the existing ensemble species distribution modelling techniques on a group of temperate bat species within various caves and mines as a proof of concept. This approach is useful in conservation practices of cavedwelling bats or any other subterranean biota [48], especially for show cave managers, as it highlights areas where anthropogenic impacts need to be minimised. The approach also offers a good insight regarding the ecology of cave-dwelling bats.

Ecology of the focal bat species
During the cold season, food scarcity may drive temperate bats (especially cave-dwellers) into prolonged torpor [77]. The optimal torpor temperatures and torpor bout durations vary according to each species' biology and fitness [27], as bats choose roosting patches that can respond to exterior temperature changes and offer stimuli if the conditions are favourable for feeding, while keeping an optimal torpor temperature in unsuitable feeding conditions [27]. If temperatures decrease, the animals become aroused and change torpor locations [21]. For example, our analysis confirmed this by the temperature range and the maximum temperature variables, which had the highest model contribution towards the mid-end hibernation period, limiting optimal habitat patches within hotter and more stable temperatures for thermophile bat species, such as Rhinolophus spp. Roost length was also relevant for bat abundance and species diversity in all of our studied subterranean environments, as Torquetti et al. [19] also found, and was linked to larger patches of stable temperatures, therefore a wider selection of optimal habitats. At the beginning of the hibernation period most species and populations were corelated with smaller distances from the entrances, because they use the climatic variability to exit daily torpor and feed in optimal conditions [78], but during the colder winter

PLOS ONE
months these correlations were weaker suggesting that populations dispersed in various optimal habitat patches, according to their body fitness [23]. Nevertheless, a general movement towards the deeper sections of the cave was observed for species such as Miniopterus schreibersii, Rhinolophus ferrumequinum or R. euryale. These populations also decreased their roosting height in the mid-end hibernation period, which suggested that the animals were searching for colder temperatures to conserve energy [79]. During the end of the hibernation season, the animals increased their activity within the roosts and move closer to the entrances, starting to feed and disperse [80]. Because of their limited fat reserves, R. hipposideros individuals preferred lower hibernation heights, where colder air masses accumulated, also observed by Ruf & Geiser [23]. The species was also described as thermophile in caves located in its northern range [78], most likely due to a lower roost stable temperature. Torpid individuals were also much more dispersed compared to other studied Rhinolophus spp. Populations during hibernation, as a high degree of clustering can increase body temperature and cause arousals [23]. This makes them much more vulnerable to anthropogenic pressures, which can be more intensely applied in the lower part of a show cave vertical section.
On the opposite spectrum of these preferences, crevice-dwelling bats such as Nyctalus noctula, or Pipistrellus pipistrellus and some populations of Myotis myotis / blythii were not pushed within the deeper parts of the roosts during the coldest periods, with some examples where populations increased their proximity to the entrances. These species have shorter torpor bout durations compared to most cave-dwelling bat species and can optimally feed in low temperatures, thus proximity to the entrances is relevant in most conditions, so the animals can naturally be aroused in optimal feeding conditions [78].
The maternity period usually requires proximity to the entrances, where temperatures are more suitable for nursing [19], except for the case of temperature inversions, where the colonies can be formed even in deep sectors of the roosts (SE 10), as our models also predicted. Often smaller sites with exterior hot air influxes are preferred [19]. Roost sectors with higher annual temperature variations were more suitable for maternity colonies (SE 3) because hot air can warm the site, offering better conditions for nursing [81]. Warmer conditions during maternity can influence embryogenesis, keeping the females more active, but also ensuring warmer nursing conditions for pups during the early stages of their life [82].
Differences in the number of individuals for each subterranean environment and activity period suggest that some roosts are used as hibernaculum (especially the larger sites-SEs 5,10) while others as maternity (smaller sites-SEs 2 and 3), as bats perform seasonal migrations between these types of roosts [83].

Species distribution models
Our approach has shown that modelling the habitat suitability of bat populations in subterranean environments is most likely less restrictive than previously considered [46], especially in large systems with more stable microclimates. External climatic factors are relevant for cavedwelling bats, influencing their movements within the roosts [16], yet we have shown that they can be easely integrated within the environmental variables. Moreover, our models were able to predict seasonal level usage of the focal species within the studied environments, due to the accurate temperature interpolations coupled with a relatively narrow niche of optimal microclimatic conditions.
The focal bat species used open spaces within their roosts to hibernate or to raise their young, allowing for a standardised occurrence sampling strategy. Data collection within a controlled modelling environment, such as a cave or a mine, ensures a consistent sampling effort, but only when dealing with open space dwellers [84]. Most troglobiontes, which make use of an extensive network of cracks and superficial habitats, would be more difficult to sample using this approach [46]. The method ensures a uniform sampling effort for most of the prediction area, offering a much better understanding of the results compared to regional or continental scale models [85].
While a smaller prediction area offers more control over the sampling strategy and the validation process, spatial autocorrelation may appear, especially when dealing with gregarious species which form colonies [79]. This can strongly influence model performance [86,87] by introducing errors resembling oversampling. However, increasing the number of surveyed sites (hence increasing the variation in the background dataset from which the models can be trained) and performing the models in multiple subterranean environments simultaneously, helped us reduce this effect, while generating valid results. Applying occurrence rarefiation methods or spatial gridded sampling strategies on our dataset, as suggested by Mammola et al. [46], to account for spatial autocorrelation would have eliminated most of the collected occurrences, generating non-valid models, but this can be further explored in larger research efforts, and especially in large subterranean environments. Heavely clustered occurrences (e.g. data from a single small habitat patch) will, nevertheless, produce extremely restrictive models, which need to be validated before continuing the analysis.
Temperature was the main limiting factor for bat colonies in subterranean environments [16], thus, our models were solely based on this parameter and had optimal validated results. This modelling approach can be more difficult in large biomass caves, especially within the tropical or equatorial regions, because the presence of bat colonies creates hotter microclimates, influencing both the environmental variables and the model results, as Lundberg & Mcfarlane showed in a study focused on temperature disturbance generated by bats [44].
Bat ranges can have a high latitude variation resulting in geographical isolated populations which are exposed to different climatic conditions [88]. For example Rhinolophus hipposideros stretches from Northern Africa to the southern part of Great Britain. This will, in turn, influence their behaviour and climatic preferences in subterranean environments, with southern populations becoming more active during the cold season [88]. Thus, our method of modelling a single taxon over a wide range of sites can lead to strong overpredictions in studies with high roost latitude variations, and needs to be further adapted and tested to account for these conditions.
Although we have used an equal contribution sum for all of the modeled species within the final SDM results, considering they are euqally susceptible to human disturbance, a further development of this approach can include weighted overlays, in which species that are found to be more susceptible to anthropogenic pressure can be prioritised in the final sensibility maps.

Implications for the conservation management of bats in subterranean environments
The models have showed abundances of favourable patches for bats in relation to their environments. We further compared these results with the anthropogenic impact to identify where conservation measures need to be applied.
Long-term anthropogenic impact can lead to roost abandonment [35], as shown in the Polovragi cave (SE 6). The cave hosted mixed maternity and hibernation colonies, but their numbers have dropped in recent years and the maternity colony has disappeared, most likely because of the upward-oriented lighting system and the human-bat interactions encouraged by the managers. The models showed sensible areas throughout the mid to end (deep) tourist sectors. Thus, a more controlled touristic access in Polovragi (SE 6) cave is required to restore bat populations, with a new lighting system pointing downward and with a strict no humanbat interaction policy. The cave should not be visited during the hibernation season and restrictions should be imposed during maternity to increase the chances of bat recolonisation.
Large bat populations can still reside in sites with a high number of tourists [78], such as the Muierilor cave (SE 5). A stable increase in temperature generated by tourists within Muierilor show cave (SE 5) was observed for the entire cold period [65] with cave lights representing a minor component. A large part of the R. ferrumequinum colony hibernated in the touristic sector. It may suggest that the animals preferred higher temperatures or were constrained to choose those areas by other unaccounted restrictions [89]. Human disturbances were strong, with halogen lighting systems projected on the hibernating colonies, noise, and flash photography. During the end of the hibernation, the colonies moved further away from the touristic sector to colder areas, most likely to conserve energy. Correlations between tourist passes and bat abundances per season were mainly not significant, suggesting that the touristic activity may be partially tolerated, as Dragu also mentioned [90]. One positive correlation was found for R. hipposideros during the beginning of hibernation; the species preferred higher torpor temperatures and maintained close proximity to the entrances, where tourists accessed the cave. This was also found by Zukal et al. who mentioned that the animals were not negatively impacted by tourism [78]. The current SDM models predicted these preferences, concentrating favourable patches in the touristic sector. Although temporary observations show no significant impact, these microclimatic preference changes can lead to an increase of artificial arousals and can on overall decrease bat body fitness on a long-term basis [79]. Specific microclimatic reconstruction projects need to be undertaken in some sectors to reduce the animal's dependency on habitat patches affected by the anthropogenic temperature changes [14]. A third small entrance in this cave that leads to a chamber that was used by hibernating bats (the Altar Hall), as proved by the guano deposits and prior observations, was closed in the past. The entrance was recently opened and then filled with rubble, but the air circulation most likely continues to flow, as ice formations can be seen near the entrance of this chamber. Sealing this opening might increase air temperature in the Altar Hall and restore a suitable area for hibernation. Bats hibernating in the touristic area might return to the previously favourable sector, decreasing the existing impact.
Valea Cetăţii show cave (SE 4) is used as a hibernaculum and harbours small bat populations. It was repurpused for tourism in 2010, and no significant population changes have been observed compared to the initial state [91]. Here, we recommend limiting cultural activities, such as concerts in the cave, during the hibernation period.
Despite the fact that Meziad show cave (SE 10) has a much larger volume than the previously described sites and it is fitted with an adequate low light system pointing downward, with no significant anthropogenic temperature spikes recorded during the monitoring period, the maternity colony from the upper level shows large population fluctuations which the touristic activity may cause. Lights in the upper level, close to the colony, should be further dimmed and managers should bypass the touristic flux in this area as much as possible to reduce stress.
The non-touristic caves located in remote regions (Stogu and Lacul Verde-SEs 7 and 8) did not show any touristic risk, but the more accessible sites, such as Gura Dobrogei cave, Canaraua Fetii mine and Hagieni mine (SEs 1, 2, and 3) were subject to vandalism and need to be closed off by special gating projects, because they concentrate large bat populations, which are crucial for the steppe bioregion. The Canaraua Fetii (SE 2) and Hagieni tunnels (SE 3) harbour large colonies and are an example of critically important bat artificial roosts. The Hagieni tunnel is mainly used for nurseries, as the SDM models and field observations suggest [92].
The Cloşani cave (SE 9) is a gated research site with controlled access. Therefore, the R. euryale hibernation colony is not submitted to significant impact, although their numbers have greatly fluctuated during recent years.
Identifying temporal and spatial bat dynamics in subteranean environments via species distribution models can help minimise accidental arousals and identify areas where ecological reconstruction techniques are needed to restore the microclimatic conditions in specific roosts [15]. As stable microclimates in subteranean environemts are changing in response to exterior temperatures [21], potentially restricting suitable habitats for cave-dwelling species [48], there is an urgent need for habitat suitability evaluations which can aid future conservation efforts. Using the existing climate change scenarios (IPCC) and the fact that temperatures fluctuations are one of the most important limiting factors for most cave-dwellers, future research can create site specific scenarios of microclimatic changes to study the distribution and occupancy of these animals [85]. Our approach can be used as a framework where species distribution modeling advancements can help decision makers apply conservation measures in light of the growing anthropogenic pressures and climate change effects, enhancing management practices for subteranean environments.